
//
// Author 
// Xianguo Lu
// lu@physi.uni-heidelberg.de 
//

#ifndef ANADEDXFIT_H
#define ANADEDXFIT_H

#define EPSILON 1e-12

class AnadEdxFit
{
 public:
  AnadEdxFit(const TString type, const Double_t p, const Double_t errlow=-999, const Double_t errhigh=-999);

  void InputData(TString histname, const TH1 *hist);

  void Fit();

  Int_t GetEntries()const;

 private:
  void SetX();
  void SetHist();
  void Save();

  //==============================================
  TH1D *fHist;

  const TString fType;
  Double_t fP;
  Double_t fErrLow;
  Double_t fErrHigh;
  Double_t fX;

  Double_t fMPV;
  Double_t fErrMPV;

  //------ fitting related ------
  TString fSFunc;
  TString fSPar[10];
  XGLUtils::FFunc fFunc; //Double_t (* fFunc)(const Double_t *, const Double_t *);
  Int_t fNPar;
  Double_t fPar[10];
  Double_t fErr[10];
  Double_t fChi[10];

  Int_t fkFindMPV;
  Int_t fkFailMPV;
  Int_t fkFailFit;
};

AnadEdxFit::AnadEdxFit(const TString type, const Double_t p, const Double_t errlow, const Double_t errhigh)
: fHist(0x0), fType(type), fP(p), fErrLow(errlow), fErrHigh(errhigh)
  , fX(-999), fMPV(-999), fErrMPV(-999)
  , fSFunc("fitting function"), fFunc(0x0), fNPar(0), fkFindMPV(-999), fkFailMPV(-999), fkFailFit(-999)
{
  const Int_t nmaxp = sizeof(fPar)/sizeof(Double_t);
  for(Int_t ii=0; ii<nmaxp; ii++){
    fPar[ii] = -999;
    fErr[ii] = -999;
    fSPar[ii]=Form("p_{%d}", ii);
  }

  SetX();
}

Int_t AnadEdxFit::GetEntries() const  
{
  if(!fHist){
    printf("AnadEdxFit::GetEntires() fHist not initialized!!\n");
    exit(1);
  }

  return Int_t(fHist->GetEntries());
}

void AnadEdxFit::InputData(const TString histname, const TH1 *hist)
{
  if(!hist){
    printf("AnadEdxFit::InputData no hist!!\n"); exit(1);
  }

  if(fHist){
    printf("AnadEdxFit::InputData fHist already exist!!\n"); exit(1);
  }

  fHist = XGLUtils::ToPDF(hist, histname);

  /*
  //test-->
  for(Int_t ii=0; ii<=hist->GetNbinsX()+1; ii++){
    printf("test %d %e %e -- %e %e %e\n", ii, hist->GetBinContent(ii), hist->GetBinError(ii), fHist->GetBinContent(ii), fHist->GetBinError(ii), fHist->GetEntries());
  }
  //test<-
  */

  SetHist();
}

void AnadEdxFit::SetX()
{
  if(fP==-999)
    return;

  if(fType.Contains("X")){
    fX = fP;
    return;
  }

  Double_t mass = -999;
  if(fType.Contains("PI"))
    mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
  else if(fType.Contains("E"))
    mass = TDatabasePDG::Instance()->GetParticle("e+")->Mass();
  else{
    printf("AnadEdxFit::SetX wrong type!\n");
    exit(1);
  }

  fX = log10(fP/mass);
}

void AnadEdxFit::SetHist()
{
  TString hname(fHist->GetName()), tit;
  if(hname.Contains("PbPb"))
    tit = "PbPb";
  else if(hname.Contains("pp")){
    tit = "pp";
  }
  else if(hname.Contains("Cosmic") && hname.Contains("h1"))
    tit = "Cosmic, dE/dx";
  else if(hname.Contains("Cosmic") && hname.Contains("h2"))
    tit = "Cosmic, dE/dx+TR";
  else if(hname.Contains("IONTR"))
    tit = "Testbeam, dE/dx+TR";
  else if(hname.Contains("ION"))
    tit = "Testbeam, dE/dx";
  else if(hname.Contains("TR"))
    tit = "TR";
  else if(hname.Contains("MC"))
    tit = "MC";
  else {
    printf("wrong hname!! in  AnadEdxFit::SetHist! %s\n", hname.Data());
    exit(1);
  }

  TString ptype;
  if(fType.Contains("E"))
    ptype = "e";
  else if(fType.Contains("PI"))
    ptype = "#pi";
  else if(fType.Contains("X"))
    ptype = "X";
  else{
    printf("wrong ftype!! %s\n", fType.Data());
    exit(1);
  }

  if(fType.Contains("X")){
    /*
    if(fErrHigh==-999)
      fHist->SetTitle(Form("%s, log_{10}#beta#gamma %.2f#pm%.2f, N %.0f", tit.Data(), fX, fErrLow, fHist->GetEntries()));
    else
      fHist->SetTitle(Form("%s, log_{10}#beta#gamma %.2f+%.2f-%.2f, N %.0f", tit.Data(), fX, fErrLow, fErrHigh, fHist->GetEntries()));
    */
    if(pow(10,fX)<100)
      fHist->SetTitle(Form("%s, #beta#gamma=%.0f, N=%.0f @", tit.Data(), pow(10,fX), fHist->GetEntries()));
    else
      fHist->SetTitle(Form("%s, #beta#gamma=10^{%.1f}, N=%.0f @", tit.Data(), fX, fHist->GetEntries()));

  }
  else{
    if(fErrLow==-999)
      fHist->SetTitle(Form("%s, %s %.1f GeV/c, N=%.0f @ log_{10}#beta#gamma %.3e", tit.Data(), ptype.Data(), fP, fHist->GetEntries(), fX));
    else
      fHist->SetTitle(Form("%s, %s %.2f #pm %.2f GeV/c, N %.0f @ log_{10}#beta#gamma %.3e", tit.Data(), ptype.Data(), fP, fErrLow, fHist->GetEntries(), fX));
  }

  fHist->SetYTitle("p.d.f.");
  fHist->SetXTitle("TRD signal (a.u.)");
  fHist->SetMaximum(15*fHist->GetBinContent(fHist->GetMaximumBin()));
}

void AnadEdxFit::Fit()
{
  if(fHist->GetEntries()==0){
    printf("AnadEdxFit::Fit no stat!! %s\n", fHist->GetName());
    return;
  }
  
  //only to test plot    
  if(0){
    //dimensionless
    fSPar[0] = "A";
    fSPar[1] = "#Delta_{p}"; 
    fSPar[2] = "#xi'"; 
    fSPar[3] = "#sigma'";
    fSPar[4] = "#kappa";
    fSFunc = "(Exp #times Landau)*Gaussian";
    
    fFunc = XGLUtils::DimlessALGE;
    fNPar = 5;
  
    fPar[0] = 9.99740e-01;
    fPar[1] = 8.91088e-02;
    fPar[2] = 1.54010e-01;
    fPar[3] = 1.79114e-01;
    fPar[4] = 6.31582e-02;
    fMPV = 2.612890e+01;
    
    for(Int_t ii=0; ii<5; ii++)
      fErr[ii] = -999;

    fErrMPV = -999;
  }

  //use ALGE fit  
  if(1){
    //dimensionless
    fSPar[0] = "A";
    fSPar[1] = "#Delta_{p}"; 
    fSPar[2] = "#xi'"; 
    fSPar[3] = "#sigma'";
    fSPar[4] = "#kappa";
    fSFunc = "(Exp #times Landau)*Gaussian";
    
    fFunc = XGLUtils::DimlessALGE;
    fHist->SetTitle(Form("%s XGLUtils::DimlessALGE", fHist->GetTitle()));
  
    fNPar = 5;

    const Double_t hmax = fHist->GetBinContent(fHist->GetMaximumBin());
    
    //------------------------------------------------
    Double_t rawpar[10];
    rawpar[0] = fHist->Integral("width"); 
    rawpar[1] = fHist->GetBinCenter(fHist->GetMaximumBin());
    rawpar[2] = fHist->GetRMS(); 

    XGLUtils::ChisquareFit(fHist, XGLUtils::ALandau, 3, rawpar, fErr, fChi, hmax/10);

    rawpar[3] = rawpar[2];    
    XGLUtils::ChisquareFit(fHist, XGLUtils::ALG,   4, rawpar, fErr, fChi, hmax/100);

    //------------------------------------------------
    
    rawpar[2] /= rawpar[1];
    rawpar[3] /= rawpar[1];
    rawpar[4] = 0;

    for(Int_t ii=0; ii<fNPar; ii++)
      fPar[ii] = rawpar[ii];

    //================= use only the very top part to estimate MPV
    fkFailMPV = XGLUtils::ChisquareFit(fHist, fFunc, fNPar, fPar, fErr, fChi, hmax/8.);

    const Double_t xmin = 0;
    const Double_t xmax = fHist->GetXaxis()->GetBinUpEdge(fHist->GetXaxis()->GetLast());
    Double_t vout[10];
    fkFindMPV = XGLUtils::GetMPVFWHM(fFunc, fPar, xmin, xmax, vout);
    fMPV = vout[0];
    
    //================= final fit ================= 
    fkFailFit = XGLUtils::ChisquareFit(fHist, fFunc, fNPar, fPar, fErr, fChi);
    
    fErrMPV = fMPV*(fErr[1]/fPar[1]);
    
    if(!fkFindMPV || fkFailMPV || fkFailFit){
      printf("AnadEdxFit::Fit fail in final LSfit! %d %d %d\n", fkFindMPV, fkFailMPV, fkFailFit);
    }
  }
  
  Save();
}

void AnadEdxFit::Save()
{
  //====================================================
  TH1D *hfit = 0x0;
  TString hname(Form("%sfit",fHist->GetName()));
  gDirectory->GetObject(hname, hfit);
  delete hfit;

  hfit = (TH1D*) fHist->Clone( hname );
  for(Int_t ib=1; ib<=hfit->GetNbinsX(); ib++){
    const Double_t ix = hfit->GetBinCenter(ib);
    const Double_t iy = fFunc(&ix, fPar);
    //printf("test %f %e -- %e %e %e %e %e\n", ix, iy, fPar[0], fPar[1], fPar[2], fPar[3], fPar[4] );
    hfit->SetBinContent(ib, iy);
  }

  //====================================================
  TH1D * hpull = 0x0;
  hname=(Form("%spull",fHist->GetName()));
  gDirectory->GetObject(hname, hpull);
  delete hpull;

  hpull = (TH1D*) fHist->Clone(hname);
  for(Int_t ib=1; ib<=hpull->GetNbinsX(); ib++){
    const Double_t cont = fHist->GetBinContent(ib);
    const Double_t iy = hfit->GetBinContent(ib);
    const Double_t ey = fHist->GetBinError(ib);
    const Double_t ip = ey>EPSILON ? (iy - cont)/ey : -999;
    hpull->SetBinContent(ib, ip);
  }

  //====================================================
  TString htit(fHist->GetTitle());
  htit = htit(0, htit.First('@'));
  TLegend *lg=new TLegend(0.18, 0.17, 0.5, 0.27);
  lg->SetFillStyle(-1);//needed to be transparent!
  lg->SetName(Form("%slg",fHist->GetName()));
  lg->SetHeader(htit);
  lg->AddEntry(fHist,"Data","p");

  if(fFunc && fPar[0]>-990)
    lg->AddEntry(hfit,Form("Fit: %s", fSFunc.Data()),"l");

  //====================================================
  TCanvas *c1 = new TCanvas;
  fHist->Draw();
  gPad->SetLogy();
  hfit->Draw("same hist C");
  lg->Draw("same");
  c1->Print(Form("plot_%s.png", fHist->GetName()));
  delete c1;

  //====================================================
  TFile *fout=new TFile(Form("hist_%s.root", fHist->GetName()),"recreate");
  fHist->Write();
  hfit->Write();
  hpull->Write();
  lg->Write();
  fout->Save();
  fout->Close();

  //====================================================
  delete hfit;
  delete hpull;
  delete lg;

  //====================================================
  //====================================================
  if(!fkFindMPV || fkFailMPV || fkFailFit){
    printf("AnadEdxFit::Save %s fitfail %d %d %d\n", fHist->GetName(), fkFindMPV, fkFailMPV, fkFailFit);
    return;
  }

  TTreeSRedirector *fStream = new TTreeSRedirector(Form("par_%s.root", fHist->GetName()));

  TString tname(fHist->GetName());
  tname= tname(tname.First('_')+1, tname.First('h')-tname.First('_')-1);//"a_IONTRh207h207pdf" -> "IONTR"

  Double_t ntot = fHist->GetEntries();
  Double_t chi2dof = fChi[0]/fChi[1];
  (*fStream)<<tname<<
    "p="<<fP<<

    "x="<<fX<<
    "ex0="<<fErrLow<<
    "ex1="<<fErrHigh<<

    "ntot="<<ntot<<
    "kfailmpv="<<fkFailMPV<<
    "kfailfit="<<fkFailFit<<
    "chi2="<<fChi[0]<<
    "ndof="<<fChi[1]<<
    "chi2dof="<<chi2dof<<
    "npar="<<fNPar;

  Double_t tmppar[10];
  for(Int_t ii=0; ii<fNPar; ii++){
    TString pn(fSPar[ii]);
    pn.ToLower();
    pn.Append("_");
    pn.ReplaceAll("'","");
    pn=pn(pn.First('#')+1,pn.First('_')-pn.First('#')-1);

    tmppar[ii] = TMath::Abs(fPar[ii]);//needed! otherwise all ps have the same value as the last one!
    (*fStream)<<tname<<
      Form("%s=", pn.Data())<<tmppar[ii]<<
      Form("e%s=",pn.Data())<<fErr[ii];
  }

  (*fStream)<<tname<<
    "mpv="<<fMPV<<
    "empv="<<fErrMPV<<
    "\n";

  delete fStream;
}

#endif
